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Abstract 

When constructing an algorithm for the numerical integration of a differential equation, 
one must first convert the known ordinary differential equation (ODE), which is defined 
at a point, into an ordinary difference equation (OAE), which is defined over an inter- 
val. Asymptotic, generalized, midpoint and trapezoidal, OAE algorithms are derived for 
a non-linear first-order ODE written in the form of a linear ODE. The asymptotic for- 
ward (typically underdamped) and backward (typically overdamped) integrators bound 
these midpoint and trapezoidal integrators, which tend to cancel out unwanted numerical 
damping by averaging, in some sense, the forward and backward integrations. 

Viscoplasticity presents itself as a system of non-linear, coupled, first-order ODEs that are 
mathematically stiff, and therefore difficult to numerically integrate. They are an excellent 
application for the asymptotic integrators. Considering a general viscoplastic structure, it 
is demonstrated that one can either integrate the viscoplastic stresses or their associated 
eigenstrains. 


INTRODUCTION 

This paper is composed of two parts. The first one discusses what we call asymptotic 
integrators, while the second one discusses their application to viscoplasticity. The purpose 
of this paper is to present new theoretical findings pertaining to these two topics. Numerical 
demonstrations are left for future papers. 
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Four, first-order, asymptotic integrators are presented, two of which are new. They are re- 
ferred to as first-order because only the linear terms are retained in the series expansions used 
in their derivations. Derivations of the asymptotic forward (7) and backward (8) integrators are 
given in Walker and Freed (1991). From these two integrators, we obtain an asymptotic, gen- 
eralized midpoint, integrator (9) — or more simply, the asymptotic midpoint integrator — which 
contains the asymptotic forward and backward integrators as limiting cases. The fourth in- 
tegrator to be considered is an asymptotic, generalized trapezoidal, integrator (10) — or more 
simply, the asymptotic trapezoidal integrator — whose derivation is given in the appendix. Like 
the asymptotic midpoint integrator, the asymptotic trapezoidal integrator reduces to the asymp- 
totic forward and backward integrators as limiting cases. After introducing these integrators, 
their asymptotic properties are discussed, and a comparison is given between them and their 
classical counterparts, viz. the forward and backward Euler, generalized midpoint, and general- 
ized trapezoidal, integrators. The first part of this paper closes with a discussion on how to apply 
the method of Newton-Raphson iteration to update the solution for the implicit integrators. 

Their are numerous viscoplastic models in the literature (see Freed et al. , 1991, for a repre- 
sentative listing) of which the vast majority are subsets of the more general structure considered 
herein. The viscoplastic theory presented below is for initially isotropic materials, such as poly- 
crystalline metals and alloys. Three types of internal stresses are considered ( i.e . the back stress, 
the drag strength, and the yield strength) of which some or all may exist in any given viscoplas- 
tic model. Their constitutive equations, which are similar in structure to the deviatoric part of 
Hooke’s law for an elastic/plastic material, are different in structure but equivalent in function 
to what one finds in the literature. This difference resides in the introduction of eigenstrains (the 
recovery mechanisms), they being governed by evolution equations, which in turn are functions 
of the viscoplastic stresses (the applied and internal stresses). Hence, one can integrate either 
the viscoplastic stresses or their eigenstrains, as we demonstrate. Which approach one ought to 
use in practice remains to be determined. 

ASYMPTOTIC INTEGRATORS 

Many physical processes, e.g. material evolution, are represented by an initial value problem 
(IVP) of a first-order ODE, as described by 

X = F[X[t}] given that X[0] = X 0 , (1) 

where X is the independent variable whose initial value at time t = 0 is X 0 . This variable may 
be scalar, vector or tensor valued in applications. The dot ‘ ' ’ is used to denote differentiation 
with respect to time. We choose time to be the dependent variable for illustrative purposes, as it 
so often is in physical applications; however, this is not a restriction on the method of numerical 
integration presented herein. We consider a particular subset of F such that 

F[X[i]] = C (*[(]] (.4 [X[t]] -X[(]) , (2) 

where the time constant C (a scalar herein) and the asymptote A (of the same type/rank as 
X) are continuous and differentiable functions of the variable X. The asymptotic integrators 
are ideally suited for equations where C > 0 and A is monotonic for a monotonically varying 
independent variable X, as is the case in viscoplasticity. If neither C nor A depends on X, 
then the equation is said to be linear; otherwise it is non-linear, as is the case in viscoplasticity. 
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Throughout this paper, square brackets [•] are used to denote ‘function of’, and are therefore kept 
logically separate from parentheses (•) and curly brackets {•} which are used for mathematical 
groupings. 

The ODE in (1) is a point function in time. Numerical integration algorithms, however, are 
based on functions evaluated over an interval in time. We therefore introduce the IVP for the 
first-order OAE, i.e. 

(3) 

X n=0 = *o , (4) 

resulting in the one-step method 


AA 


= T [X nU ] 


where 


AX = X n+ i - X n with 


X n+ i = X n + 7- [X n+ ^]- At . 


(5) 


This relationship is used to represent the IVP of the first-order ODE in numerical integration. 
We use the notation X n = X[t], X n+( ^ = X[t+<pAt] and X n +i = X[t+At], where 0 < 0 < 1. 
The formulation is explicit whenever 0 = 0; whereas, it is implicit whenever 0 < 0 < 1. Our 
objective is to obtain 7 for the forward, backward, midpoint, and trapezoidal, asymptotic 
integrators. 

One can introduce an integrating factor into the differential equation (1—2), and as a conse- 
quence obtain the recursive integral equation 


[ rt+At 

-J c[xi(\)dt x[t] 

+ jf' + “exp[-jf‘_ + "c[X[C]l<iC C[X[fMX[£M 


( 6 ) 


which is the exact solution to this first-order ODE. John Bernoulli (1697) developed a non- 
recursive solution similar to (6) for the equation X = &X + bX n , whose solution was sought 
earlier by Jacob Bernoulli (1695). John’s solution is expressed as a quadrature, since the integral 
of dzjz in the form of a logarithm was not generally known until later that same year (cf. Ince, 
1956). The second line in (6) accounts for the non-homogeneous contribution to the solution. 
It is a Laplace (1820) integral when the integrand has its largest value at the upper limit t+At, 
and therefore possesses an evanescent memory of the forcing function C[X[£]].A[X[£]] provided 
that C[X[C]] > 0 over the interval (M + At). This fading memory means that the solution 
will depend mainly on the recent values of the forcing function, and that by concentrating the 
accuracy on the recent past we obtain accurate asymptotic representations of the solution. 

Walker and Freed (1991) obtained a variety of asymptotic and exact solutions to (6) by 
expanding the functions C and A into series, and then integrating term by term; in particular, 
one obtains the first-order, asymptotic, forward integrator, 


X n+1 = e~ c ^- At X n + (l - e - c l x " 1At ) A\X n \ , 


(7) 


and the first-order, asymptotic, backward integrator, 

X n+1 = e-W+il-At Xn + (l - e- c [ x "+d-At) A { Xn+1 } , (8) 
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by expanding C and A in Taylor series about the lower and upper limits of integration, respec- 
tively, and truncating them after their linear terms. Application of these integrators to perfect 
plasticity has been discussed by Simo and Taylor (1986), and used by Margolin and Flower 
(1991) who found asymptotic integration to be more accurate than the radial return algorithms. 
Application has also been made to viscoplasticity by Walker (1981, 1987), Ju (1990), Walker 
and Freed (1991), Iskovitz et al. (1991), and Freed and Walker (1992), who have shown that 
asymptotic integration is preferred over classical integration for viscoplastic analysis. 

An asymptotic midpoint integrator (a new algorithm) is easily acquired by adding the con- 
tribution from the implicit integrator (8) taken over the interval (t,t + <pAt) with that of the 
explicit integrator (7) taken over the remaining interval (t+0At, t-f At), resulting in 

X n+1 = e~ c ^^ At X n + (l - e- c [* B +*]- A< ) A{X n+ t) . (9) 


Whenever 0 = 0, this integrator reduces to the forward integrator (7); likewise, whenever 0 = 1, 
it reduces to the backward integrator (8). 

The asymptotic trapezoidal integrator (another new algorithm) is obtained from the same 
procedure used to acquire the forward and backward integrators, except now the series expan- 
sions are weighted about both limits of integration, instead of just one. The resulting integrator 
is 


X„ + . = e -* 1 *"*-!' 41 X. + (l - e-W" A t [X„, X n+1 ] , 


where 


( 10 ) 


A 1 [X n ,X n+1 ] = (1 -<f>)A[X n } + <t>A[X n+1 ], (11) 

Cx[X n , X n+l ) = (1 - 0)C[X m ) + 0C'[X w+1 ] , (12) 

C 2 [X n ,X n+l ] = (1 — 4>) 2 C[X n ] + 0(2 — <p)C\X n+ ]\ , (13) 

and whose derivation is given in the appendix. This integrator distinguishes itself from the prior 
three integrators in that it is a two-point integrator, i.e. it has a functional dependence on both 
X n and X n +i, instead of a one-point integrator. Whenever 0 = 0, this integrator reduces to 
the forward integrator (7); likewise, whenever 0 = 1, it reduces to the backward integrator (8). 

These four asymptotic integrators are exact and equivalent whenever C and A are both 
constants, i.e. whenever the ODE is linear; otherwise, they are approximations. Higher-order 
asymptotic integrators can be found in Walker and Freed (1991). 

Asymptotic Properties 


To acquire the asymptotic backward integrator, the integrands in (6) were expanded in 
Taylor series about their upper limits (Walker, 1987, and Walker and Freed, 1991), where each 
integrand has its largest value and contributes the most to the integral. By retaining but a 
single term in the Taylor series expansions, the integrands are accurately approximated where 
they are largest, and the neglect of the higher-order terms is only felt near the lower limits where 
each integrand contributes only a small amount to the integral because of its exponential decay 
from the upper limit. The neglect of the higher-order terms in the Taylor series thus results 
in an algorithm that is asymptotically correct at the upper limit; in particular, the asymptotic 
expansion for the backward OAE (8) is given by 

lim X n+i x A[AT n+ i] , (14) 

At— +large 
C> 0 
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which is the correct asymptote for the ODE in (1-2). Normally, when treating asymptotic 
expansions, the exponential decay of the integrand allows the lower limit to be replaced with 
zero or minus infinity to ease the integration. This was not done in the present case, however, 
so that by retaining the lower limit as t, we obtain a uniformly valid asymptotic algorithm in 
the fully implicit approximation, provided that C[X[(]] > 0 over the interval (t,t-\-At). 

To acquire the asymptotic forward integrator, the integrands in (6) were expanded in Taylor 
series about their lower limits (Walker and Freed, 1991), where the neglect of the higher-order 
terms in the Taylor series results in integrands that become progressively more inaccurate as 
they approach their upper limits where the contribution from each integrand is most important. 
The explicit approximation is not, therefore, a valid asymptotic representation of the integral 
when the Taylor series is truncated at a finite number of terms. In contrast with (14), the 
asymptotic expansion for the forward OAE (7) is given by 

lim X n+ i x A[X n \ , (15) 

At — ►large 
C>0 

which can produce the correct asymptote for the ODE in (1-2), but only when A[X n ] « 
A[X„ + i]; in other words, X n must be in the neighborhood of its asymptote A in order for the 
asymptotic forward integrator to predict an accurate asymptotic value. 

The asymptotic midpoint and trapezoidal integrators contain the forward and backward 
asymptotic integrators as their limits. Their respective asymptotic expansions are 

lim X n+ i x A[X n+<t> ] , (16) 

At — ►large 

C> 0 


and 


lim Jf x 

At— ►large 

< 7>0 


(Ci [Xn 1 Xn ± l\ 

\C 2 [X n ,X n+1 ] 



*»+l] 


> 


(17) 


which, like the asymptotic forward integrator, can both produce the correct asymptote for 
the ODE in (1-2), but only when A\X n+<j> ] « A[X n+ i] or when C[X n ] « C[X n+1 ] and 
A[X n 1 * A[X n+1 ], respectively. The asymptotic properties of these two integrators improves 
as 0 approaches the value of 1. Why then would one choose a value for 0 other than 1? The 
answer is because when 0 = 0, these integrators produce an underdamped transient response; 
whereas, when 0=1, they produce an overdamped transient response. Therefore by choosing 
an appropriate value for 0, which will depend on the specific ODE being integrated, one can 
effectively cancel out any unwanted damping in the solution that is caused by the integrator. 
Unfortunately, there is no known rigorous procedure for determining an optimum value for 0. 

All four asymptotic integrators have the desirable property that X goes to A as the ODE 
saturates, independent of whatever numerical error may have accumulated while integrating 
through the transient domain. This property is not shared by their classical counterparts. 
Multiple time steps are required to mitigate oscillations when integrating to saturation using 
0 < 0 < 1; remarkably, one time step is often sufficient when integrating to saturation using 
0 = 1 . 


Asymptotic vs. Classical OAEs 

Rearranging the recursive solutions of (7-10) into the form of the difference equation in (5), 
one obtains the one-step, asymptotic, forward integrator, 


X n+l = X n + F f [X n ]-A t y 


(18) 
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where 


/I _ e -C[X,]-A*\ 

^J = ( wr) Wl 


the one-step, asymptotic, backward integrator, 

X n +\ = X n + F b \X n+1 ]-At , 

where 


Fy[X n+1 ] = 


* g'\~C\X n + l] 


1 


C[X n+1 ]-At ) 
the one-step, asymptotic, midpoint integrator, 

X n +1 = X n + J- ,f\X n^- At , 


F[X n+1 ] ; 


where 


given that 


- ( ! _ „ (! _ e-cn.wl-4.) J 1 ” + * 1 


X n +^ s= X n + <pT b [X n+< t,]-At 

= x n + <p 


e +*C[X n+ 4] At _ X \ 


F[X n+ +]-At 


<t>C[X nU YAt ) 

± (1 - (f))X n + <t>x n+1 unless 0=1; 

and the one-step, asymptotic, trapezoidal integrator, 

^n+1 = X n +F t [X n ,X n+l ]-A t, 


(19) 

( 20 ) 
( 21 ) 

( 22 ) 

(23) 

(24) 

(25) 


where 


F t [X n ,X n+1 ] 


C x [X n ,X n+1 ] 


(( 


g £*2 [Xn ,JT n + l l'A< 

C 2 [X n ,X n+1 ]-At 


^jA l [X n ,X n+l ] 


^ [-X'n jXn -f 1 j ' 

C 1 [X n ,J»C n+1 ]A( 



(26) 


which, unlike the prior ^s, cannot be written as a simple function of F. It is easily verified that 
(18) and (20) are the Umiting cases of (22) and (25). All four asymptotic difference functions 
T are non-linear and have m explicit dependence on the time step size At. 

Contrast these asymptotic integrators with their classical counterparts; namely, the forward 


Euler integrator, 


X n+1 = X n +F f [X n ]-A t, 

(27) 

where 


F f [X n ] = F[X n ] ; 

(28) 

the backward Euler integrator, 


-Xn+l = X n + JF 6 [X n+ i]-At , 

(29) 
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where 


(30) 


T b \X.^\ = F\X ntl \ ; 

the classical, generalized, midpoint integrator, 

X n +\ = X n + 7F ^[Xn+j]- At , (31) 

where 

= n*.+*l ( 32 > 

given that 

Xn+j, = X n + <pFb[X n +$\ "At 
= X n + <f>F[X n+</> )-At 

= (1 -0)X n + 0X„ + i V 0; (33) 

and the classical, generalized, trapezoidal integrator, 

X n+l = X n +F t [X n ,X n+1 ]-At , (34) 


where 


r t [x n ,x n+1 } = (i - 4>)r f \x n ] + <pr„[x n+l } 

= (1 — 4>)F[X n ] + 0F[X n+1 ] . (35) 

Obviously, (27) and (29) are the limiting cases of (31) and (34). Unlike their asymptotic coun- 
terparts, the classical difference functions T are linear and do not have an explicit dependence 
on the time step size At. 

The classical midpoint integrator (31) is unconditionally stable over the interval l / 2 < 4> < U 
curiously, the classical trapezoidal integrator (34) is not stable at 0 = l / 2 for non-linear problems 
(Hughes, 1983). With respect to a local truncation error, they are second-order accurate when 
0 = i/ 2 . Similar statements pertaining to the properties of stability and truncation error, as 
they apply to the asymptotic midpoint and trapezoidal integrators (22 and 25), are not presently 
available. 

The extra terms in the asymptotic difference functions (19, 21, 23), the likes of which are not 
found in their classical counterparts (28, 30, 32), are the source of the asymptotic characteristics 
of our integrators; characteristics not shared by the classical integrators. These asymptotic 
integrators reduce to their classical counterparts only in the limit as C At — * 0. The asymptotic, 
trapezoidal, difference function (26) bears little resemblance to its classical counterpart (35). 

Midpoint Solution Procedure 

To effect the integration — via the midpoint integrators — for updating the solution from X n 
to X n+ i, one can adopt the following procedure. First, use a backward integrator (24 or 33) 
over the interval (M + 0Af) to determine X n+4> . Then, use the midpoint integrator (22 or 
31, respectively) over the interval (t,t + At) to determine X n+1 . In essence, we are implicitly 
integrating over the interval (M+0Ai), and then explicitly integrating over the interval (t+At) 
at some midpoint t+4>At. 
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Applying the method of Newton-Raphson iteration to enable the implicit integration of (24 
or 33) over the interval (t, t+<f>At) leads to a correction 6X for the z‘ h iterate approximation 
that is resolved from the linear system of equations, 

|/ - 0At ) | -ex = X n + - X { ;\, , (36) 


where I is the identity matrix. The next approximation is then calculated from 

= X« * + SX . (37) 

which is substituted back into (36) with iteration continuing until the convergence criterion, 


ll-x. + WdVAf - *i+4 < e , 


(38) 


is satisfied for some norm || • || and maximum tolerance or error e. Now that X n+tj> is a known 
quantity, one uses the associated midpoint integrator (22 or 31) over the interval (t,t+At) to 
update X n+1 . 

This algorithm, although appearing to be different, is equivalent to the one used by others 
(see, for example, Hornberger and Stamm, 1989) who, in effect, combine the implicit and mid- 
point integrations into the Newton-Raphson iteration. One can combine these two steps into a 
single one for the classical midpoint integrator, because of the linearity in (33); however, this 
cannot be done for the asymptotic midpoint integrator, because of the non-linearity in (24). 

The derivative of in (21) with respect to X n+ ^ is the Jacobian in (36) for the asymptotic 
midpoint integrator; it is, 


= f e+^lO-* - 1 1 ( 9A[X^ \ 

dX n ,<, \ J 1 ” + *'( dX„ +t ~ J 

+ A ‘} (A[Xi%] - X!H 4 ) , (39) 


which has an explicit dependence on the time step size <pA t. Similarly, from (30) one obtains 
the Jacobian associated with the classical midpoint integrator; it is, 




dX 


n+<f> 


= C[X 


(«) 1 

n+<f > . 






dX 


-I + 


dX 


n+<j> 




+4> 


X (i) ) 
-A-n+t) 


(40) 


which contrary to (39), does not explicitly depend on the time step size 4>A t. The asymptotic 
Jacobian differs from the classical one by the coefficients found in the curly brackets of (39). 
Both of these Jacobians are expressed in terms of the derivatives of the time constant C and 
the asymptote A taken with respect to the independent variable X n+< p. 

Trapezoidal Solution Procedure 

The updating of the solution from X n to X n+ i using the trapezoidal integrators (25 or 34) 
can also be accomplished through Newton-Raphson iteration, where 


I - At 




dX 


n + 1 


■6x = x.+^[x.,xf»,]-A t - x;;;, , 


'(>) 


(41) 
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( 42 ) 


with the next approximation being calculated from 

x!S;> = *2., + ex , 

which is substituted back into (41) and iterated upon until the convergence criterion, 

llXn + FtlX^X^U-At-XtU <e, (43) 


is satisfied; whence, the solution is updated. 

The derivative of Ft in (26) with respect to X n +i is the Jacobian in (41) for the asymptotic 
trapezoidal integrator; it is, 


dX n+l 


<P 


Hc 2 [X„,x2il'A</ l 8X "+> 

,,, _S4[x( i >,l\ ^ rv vO) 1 A * d( 


-b) l 
n+lj 


v \ n* f*Tli / \ 

+ £,[*„, X„ t ,| gx ^ j e 9X n+1 A„ 

/ . AT r V v(*) 1 A 4 _ \ 


• t)A„ + , j 

,J(i + c 2 [x„,x2ilAt) 

+ ( (c 2 [x„,x« 

r -..ftl i 


-CtlX^hYAt _ 
n’i “ntu / 

(C 2 [X„,X2,]A1) 2 


which has an explicit dependence on the time step size At. Similarly, from (35) one obtains the 
Jacobian associated with the classical trapezoidal integrator; it is, 


ar.lx^xS,! 

3X„,, 


= 4 



( M[x 2 il j 
1, ax„ +1 



(45) 


which contrary to (44), does not explicitly depend on the time step size At. Both of these 
Jacobians are expressed in terms of the derivatives of the time constant C and the asymptote 
A taken with respect to the independent variable X n+ i. 


VISCOPLASTICITY 


The analysis of metallic response for high temperature applications requires mathematical 
models capable of accurately predicting short-term plastic strain, long-term creep strain, and 
the interactions between them. Viscoplastic models attempt to do that. Multiaxial, cyclic 
and non-isothermal loading histories are normal service conditions — not exceptional ones — all 
of which challenge the predictive capabilities of such models. 

The stress a tj is taken to be related to the infinitesimal strain e through the constitutive 
equations of an isotropic Hookean material, viz. 

Sij = 2/x(£ 0 - - e?) where 4 = 0 , (46) 

and 

a tl = 3 K(e tt - a(T - T 0 )4 ) , (47) 
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which are characterized by the shear // and bulk k elastic moduli, and where 


Sij = ~ Vs °kk and E l} = e X} - 1/3 e kk ^ (48) 

denote the deviatoric stress and strain, respectively. The constant a is the mean coefficient of 
thermal expansion, which acts on the difference between the current temperature T and some 
reference temperature T 0 . The quantity 6 is the Kronecker delta; it is 1 if i = j, otherwise it 
is 0. Repeated Latin indices are summed over from 1 to 3 in the usual manner. Equation (46) 
characterizes the deviatoric response, while (47) characterizes the hydrostatic response. The 
plastic e? and thermal a(T — Tq)^ strains are eigenstrains that represent deviations from 
deviatoric and hydrostatic elastic behaviors, respectively. 

A typical viscoplastic model can admit any of three, phenomenological, internal, state vari- 
ables; they are: i) a deviatoric back stress B it) a drag strength D, and Hi) a yield strength 
Y. The back stress accounts for kinematic hardening effects, while the drag and yield strengths 
account for isotropic hardening effects. These internal variables are described through their 


constitutive equations 




Ba = 2 H(e%-X i} ) 

where 

= 0 , 

(49) 

D = D 0 + h(p-x) 

where 

0 

A 

0 

Q 

(50) 

Y = Y 0 + h(p-x) 

where 

T 0 >0, 

(51) 


and where H, h, and h are positive- valued hardening moduli; X ijx an( f X are eigenstrains; D 0 
and T 0 are the annealed (minimum) values for the drag and yield strengths; and 


V = 



de? 

dr 


dr , 


(52) 


is the accumulation of plastic strain incurred over the loading history. Notice the similarity 
between the constitutive equation for stress (46) and those for internal state (49-51). In Hooke’s 
law, the eigenstrain s? represents a deviation from elasticity — a phenomenon called plasticity. In 
the constitutive equations for internal state, the eigenstrains X t] , x, and x represent deviations 
from strain hardening — phenomena called recovery. This competition between the mechanisms 
of strain hardening and recovery is in accordance with Bailey’s (1926) hypothesis that there is 
a “balance between the generation and removal of strain-hardened material” at steady state. 

The eigenstrains in the constitutive equations (46, 49-51) are governed by evolution equa- 
tions. In our general viscoplastic structure, the plastic strain is considered to evolve according 
to the Prager (1949) relation, 

j|5§| , (53) 

where l / 2 {Sij - B l} )/\\S - B\\ defines the direction of plastic flow. The back stress is taken to 
recover according to the Armstrong and Frederick (1966)/Chaboche (1977) relation, 

= V 2 (“ llil + tftf) jj§| , (54) 

where - l / 2 B t j/\\B\\ defines the direction of recovery for the back stress. Here L > 0 is the 
limiting state of back stress, > 0 is typically an Arrhenius function of temperature, and 
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R > 0 is the thermal recovery function. The drag and yield strengths are considered to recover 
according to the Orowan (1947)/Ponter and Leckie (1976)/Chaboche (1977) relations, 


X = f II i'll+tfr, (55) 

& = jll e p || + 1?r , (56) 


where £ > 0 and t > 0 are the limit functions, and r > 0 and f > 0 the thermal recovery 
functions. In these evolution equations for the eigenstrains of internal state (54-56), the first 
term on the right hand side accounts for dynamic or strain-induced recovery, while the second 
one accounts for thermal or time-induced recovery. 

The kinetic equation of evolution is taken to be governed by a Zener and Hollomon (1944) 
decomposition in state; that is, 


||e p || = tf[T] Z 


|S-B|| -r' 
D 


(57) 


where the Zener parameter Z > 0 is a temperature normalized magnitude for the plastic strain 
rate, which we assume incorporates a Bingham (1916)/01droyd (1947)/Prager (1949) criterion 
for yield. The Macauley bracket (x) has the value of 0 whenever x < 0, or the value of x 
whenever x > 0. 

The norms (or magnitudes) pertaining to the deviatoric tensors used herein are defined by 
\\I\\ = y/Vilij ~hj and || J|| = T X] J tj , (58) 


where I tj is any deviatoric stress-like tensor, and J i} is any deviatoric strain-like tensor. These 
are the norms of von Mises (1913), where the coefficients have been chosen to scale the theory 
for shear. 

The constitutive, evolution, and kinetic equations (49-57) become specific, thereby defining 
a particular viscoplastic model, when material functions for Z , H , h, h , L , £, £, R, r, and f 
are provided. 

Asymptotic Form for Eigenstrains 

To gain a solution for this viscoplastic theory, one may numerically integrate the evolu- 
tion equations for the four eigenstrains. This is accomplished by substituting the constitutive 
equations (46, 49-51) into their respective evolution equations (53-56), thereby obtaining 


whose time constants are given by 


4 

= C”(Al- 

4) > 

(59) 


l 

* 

o 

ii 

- X v) . 

(60) 

X 

= c x (a x - 

x) , 

(61) 


1 

•X 

-IS 

'X 

o 

II 

x) , 

(62) 

— 

»+ H 
S-B 11 1 

1 >0, 

(63) 
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( 64 ) 


C * = ~ e \m + ^ r >0, (65) 

C* = )l|i p ll + £^><>, (66) 

and whose asymptotes are given by 

A -> - iir (67) 

4-5 = 4 . ( 68 ) 

A* = ^2+p, (69) 

A* = 5+j>. (70) 

h 


These are well behaved equations for the asymptotic integrators because i) the time constants 
are non-negative, and ii) the asymptotes are monotonic and bounded for monotonically varying 
independent variables. The inequalities associated with the time constants require that ||e p || = 0 
when — B\\ = 0, that R = 0 when ||S|| = 0, and that f = 0 when Y = Yq = 0, all of which 
are reasonable restrictions. 

Asymptotic Form for Stresses 

An alternative means of updating the solution is to integrate the viscoplastic stresses rather 
than their eigenstrains. This is accomplished by differentiating the constitutive equations 
(46, 49-51) and then combining these relations with their evolution equations (53-56), thereby 
resulting in the more complex system of equations 


Si j = C s (4 - s.,) , 

(71) 

B,j = C B (A* - B it ) , 

(72) 

b = c d (a d -d ) , 

(73) 

Y = C Y (A r - Y) , 

(74) 


whose time constants are given by 


C s 

C B 

C D 

C Y 




H 


h 


S-B 

' L + \\S-B 




|e p 


L\\S-B\\ 

fir 


+ 


D 


[ ) Hi' 

h 2 ) ’ 


+ 


tiR 

\\B\\ 




(75) 

(76) 

(77) 

(78) 
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and whose asymptotes are given by 


,s = 2E„ + (llel/115-Bll)g 0 - 

Ums-BW-m/dry^t' 

, fl _ (||e p ||/||S-B||)^ 

((L + ||5-B||)/L||5-B||)||el -H>*/||B|| - H/F 2 ’ 

aD = ii* p ii - (vm 

H^ll /Utir/D-h/h 2 ’ 

A v = ll£ p ll-(A/h 2 )r 0 

||e p ||/^ + tff/T- A/h 2 ’ 


(79) 

(80) 
(81) 
(82) 


These relationships account for the fact that the shear modulus varies with temperature, and 
that the hardening moduli H, h , and h may vary with state, as they do in some viscoplastic mod- 
els; hence, their rates of change must be taken into consideration when integrating the stresses. 
This complexity is not present in the previous approach where the eigenstrains are integrated. 
Because of these rates of change in the moduli, it is possible — although not probable that 
the time constants C may become negative valued, and as a consequence, the asymptotes may 
become unbounded. (Unlike the previous approach, here the time constants appearing in the 
denominators do not divide out in the equations for the asymptotes; hence, the capability of 
these asymptotes becoming unbounded.) Nevertheless, integrating the stresses does have one 
advantage over integrating the eigenstrains; that is, the asymptotes Afj, Af jt A D , and A Y be- 
come stationary at steady state; whereas, A p {j , A*, A x , and A x continue to vary at steady state. 
In all of our implementations of asymptotic integrators applied to viscoplastic models (Walker, 
1981, 1987, Walker and Freed, 1991, and Freed and Walker, 1992), the viscoplastic stresses were 
always the integrated variables. 

SUMMARY 


This paper presents two, new, numerical integrators — the asymptotic midpoint and trape- 
zoidal integrators — which contain the asymptotic forward (0 = 0) and backward (0 = 1) inte- 
grators as their limiting cases. By varying 0 between 0 and 1, one can minimize any unwanted 
damping in the solution that is caused by the integrator. The forward integrator is typically 
underdamped, while the backward integrator is typically overdamped. Which of these two, new, 
asymptotic integrators is preferred for applications is a subject of future work. 

Viscoplasticity is an application where the asymptotic backward integrator has been used 
with great success. A general viscoplastic structure is presented in this paper. Into this structure 
we introduce a new concept referred to as the viscoplastic eigenstrains, which account for various 
recovery /softening mechanisms. This concept enables one to better visualize the structure of 
viscoplastic theory, with an added benefit of providing two options for updating/integrating the 
viscoplastic IVP. One may either integrate the eigenstrains or their viscoplastic stresses. Both 
approaches have their plus and minus points from a theoretical perspective, but whether one 
ought to integrate the eigenstrains or their stresses in practice is a subject of future work. 
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APPENDIX 


The derivation of the asymptotic trapezoidal integrator (10) is given below. The method of 
solution is: i ) obtain an appropriate series expansion for the integral in the two exponentials 
of (6), ii) expand j 4[X[£]] and C[X[f]j into appropriate series expansions, and Hi) analytically 
integrate the remaining integral— the non-homogeneous contribution— with its approximated 
integrand. The trapezoidal algorithm, by definition, weights the integrand at its two limits of 
integration. 

For the first step in the solution procedure, let us consider a Taylor series expansion for an 
integral about it lower limit, 


[ /[£] = f[ a \{b - a) + y 2 /[a](& - a) 2 + 1 / 6 /[a](6 - a) 3 

Jf=a 


+ 


(83) 


)(,=a 

and one about its upper limit, 


/ 6 m d£ = f[b](b - a) - y 2 f[b](b - af + i/ J[b)(b - a) 3 , (84) 

which assumes the integrand / to be continuous and differentiable over the domain (a, h). Trun- 
cating these expansions after their linear term, taking 1-0 of the lower one and 0 of the upper 
one, where 0 < 0 < 1, and then adding these two contributions, allows the integral in the two 
exponentials of (6) to be approximated as 

X , ~ _-{(W)C[x n ]+*c[x n+ i]}-At x 

+ [ t+ *\-{V-M[m+MXn+m+*-Q C[X[£]] A[X[f]]d$ , (85) 

Jz=t 
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given that the time constant C is continuous and differentiable over the time domain (t,t+At). 

In the second step of the solution procedure, the functions .A[X[£]] and C[X[^]] are expanded 
in Taylor series about t and t+At, which bound £, and are then averaged such that 

X[XK]] = (i-0){4[x„] + A[x„]«-i) + iM[x n ]«-() 2 + ---} 

+ 0{^[x„ + ,]-A|jr„ + ,](( + Ai-o + iM[x„ tl ](( + A(-fl 2 --- } ,(86) 

and 

C[XK]| = (1 - 0) (qx n ) + C|X„]($ - () + i/ 2 c[x„]K-() 2 + ---} 

+ <j> {c[X — C[X n+ i](t + A1 — £) + + A( — £) 2 — • • ■} . (87) 

Truncating each expansion after its linear term enables the integrand of (85) to be further 
approximated as 


X , , ~ e -{(l-<t>)C\X„\-\-4>C\X nJri ])-^t jr 

Ti+l — ° ^ n 

+ ((1 - <f>)C[X n ] + <pC{X n+1 }) ((1 - <P)A[X n ) + 4>A[X n+1 ]) 

rt+At 

X / e -{(l-^C[X n ]+0(2^)C[X n+ i]}( t+ A t -O ^ 

h=t 


( 88 ) 


This step requires both the asymptote A and the time constant C to be continuous and differ- 
entiable over the time interval (t,t+At). 

The final step is to analytically integrate the remaining integral in (88). To do this, we 
introduce the change in variable Z — t + At — £, and hence rewrite (88) as 


X n+1 ~ e -{(l-^)C[X„l+^C[X„ + 1 ]>A* Xw 

+ ((1 - <f>)C[X n ] + 4>C[X n+l ]) ((1 - <j>)A[X n ) + <f>A[X n+1 ]) 

x [** g-{(l-^ 2 C[X„]+^(2-^)C[X n+1 ]}Z dz 

Jz = 0 

such that upon integration one obtains 

X n + 1 ~ g-{(l-^)C[X»]+^C[X n+ i]} A« Xn + ^ _ e -{(l-^) 2 C[X n ]+^(2-^)C[X„ + 1 ]}.Atj 

V(l-<*) 2 C[X„]+0(2-0)C[X„ + ,]J“ («^[X„]+M[X„ + ,1) 


which is the asymptotic trapezoidal integrator given in (10). 


(89) 


(90) 
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